clear all;
close all;
delete(gcp('nocreate'));
%tic
global usual_delta;
usual_delta = 0.05;
% Current folder for Matlab codes
fp.matlab = pwd;
common_path = extractBefore(fp.matlab,"\matlab_dir");

% Paper folder;
fp.paper = sprintf('%s\\paper', common_path);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%fixed parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%We estimate the model at lambda=0.95;
fp.lambda = 0.95;

%model primitives;
fp.ages =9:12;
fp.periods=length(fp.ages);
fp.percentiles_kid = linspace(5,95,10);
fp.hc_points_kid = length(fp.percentiles_kid);

fp.Min_Time = 0.01;
fp.Max_Time= 1;
fp.inv_points=20;
fp.percentiles_inv= linspace(5,95,fp.inv_points);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%indexes for latent data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fp.ind_data.child_id = 1;
fp.ind_data.child_age = 2;
fp.ind_data.n_sim = 3;
fp.ind_data.child_C = 4;
fp.ind_data.child_C_tp1 = 5;
fp.ind_data.inv = 6;
fp.ind_data.peers = 7;
fp.ind_data.peers_log = 8;
fp.ind_data.ParStyle = 9;
fp.ind_data.peers_tp1 = 10;
fp.ind_data.peers_log_tp1 = 11;
fp.ind_data.neighborhood = 12;
fp.ind_data.moved = 13;

fp.ind_data_network.child_id=1;
fp.ind_data_network.child_age=2;
fp.ind_data_network.n_sim=3;
fp.ind_data_network.child_C=4;
fp.ind_data_network.n_friends=5;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Loading Estimated Moments;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Loading Initial Conditions (Table 4);
initial_conditions = importdata('initial_conditions.txt');    
fp.schools_distrib = initial_conditions;

%Loading Income information of neighborhoods (it affects preferences for parenting in the model);
mean_income_by_neighborhoods = importdata('family_income_neighborhoods.txt');
mean_income_by_neighborhoods = mean_income_by_neighborhoods./10000;
fp.mean_income_by_neighborhoods = mean_income_by_neighborhoods;


%loading estimated elasticities of neighborhood wrt school quality (from
%Agostinelli, Martellini and Luflade, 2022)
elasticity_moving_by_kid_skills=importdata('elasticities_by_kid.txt');
fp.elasticity_moving_by_kid_skills = elasticity_moving_by_kid_skills;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%random mumber draws
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

rng(10); %fix seed;
fp.n_rip_sim = 200;                         % # of simulation of the economy;
fp.tot_draws_mcintegration_network = 50 ;   % # of simulation of the network;

%Initial skills;
fp.init_draws{1} = norminv(rand(fp.schools_distrib(1,3),1,fp.n_rip_sim));
fp.init_draws{2} = norminv(rand(fp.schools_distrib(2,3),1,fp.n_rip_sim));
fp.init_draws{3} = norminv(rand(fp.schools_distrib(3,3),1,fp.n_rip_sim));
fp.init_draws{4} = norminv(rand(fp.schools_distrib(4,3)+50,1,fp.n_rip_sim));

%Peer shocks (backward induction)
fp.peers_backward{1} = rand(length(fp.percentiles_kid)*length(fp.percentiles_kid)*fp.inv_points,fp.schools_distrib(1,3),fp.tot_draws_mcintegration_network, fp.periods-1);
fp.peers_backward{2} = rand(length(fp.percentiles_kid)*length(fp.percentiles_kid)*fp.inv_points,fp.schools_distrib(2,3),fp.tot_draws_mcintegration_network, fp.periods-1);
fp.peers_backward{3} = rand(length(fp.percentiles_kid)*length(fp.percentiles_kid)*fp.inv_points,fp.schools_distrib(3,3),fp.tot_draws_mcintegration_network, fp.periods-1);
fp.peers_backward{4} = rand(length(fp.percentiles_kid)*length(fp.percentiles_kid)*fp.inv_points,fp.schools_distrib(4,3)+50,fp.tot_draws_mcintegration_network, fp.periods-1);

%Peer shocks (forward induction)
fp.peers_forward{1} = rand(fp.schools_distrib(1,3),fp.schools_distrib(1,3),fp.n_rip_sim,fp.periods-1);
fp.peers_forward{2} = rand(fp.schools_distrib(2,3),fp.schools_distrib(2,3),fp.n_rip_sim,fp.periods-1);
fp.peers_forward{3} = rand(fp.schools_distrib(3,3),fp.schools_distrib(3,3),fp.n_rip_sim,fp.periods-1);
fp.peers_forward{4} = rand(fp.schools_distrib(4,3)+50,fp.schools_distrib(4,3)+50,fp.n_rip_sim,fp.periods-1);

%Parenting Style shocks
fp.PS_forward{1} = rand(fp.schools_distrib(1,3),fp.n_rip_sim,fp.periods-1);
fp.PS_forward{2} = rand(fp.schools_distrib(2,3),fp.n_rip_sim,fp.periods-1);
fp.PS_forward{3} = rand(fp.schools_distrib(3,3),fp.n_rip_sim,fp.periods-1);
fp.PS_forward{4} = rand(fp.schools_distrib(4,3)+50,fp.n_rip_sim,fp.periods-1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%COUNTERFACTUALS;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fp.origin_neigh         = 1;
fp.receiving_neigh      = 4;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Simulate Baseline Economy        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

load('param_est');   
fp.do_estimation =0;
fp.do_counter = 0;            
fp.do_counter_type = 0 ;

obj_eval = obj_function( param_est, fp );
data_baseline_origin = obj_eval.simul_data{fp.origin_neigh};

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%Counterfactuals;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;

%Here we double the shocks for peers formation in the backward problem
%because in the counterfactuals we double the state space (moving families
%with different preferences into the new neighborhood);

fp.peers_backward_original = fp.peers_backward;

for j = 1:1:4
   
    fp.peers_backward{j} = [ fp.peers_backward{j} ; fp.peers_backward{j} ] ; 
 
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Policy #2: Policy Effects of moving average many children from Sending Neighborhood)
%           allowing for parents in receiving neighborhood to leave;

%Appendix Figure E-4: Residential Elasticities w.r.t. Mean School Peer Achievement
%Appendix Figure E-5: Policy and Scaling Effects on Skills (with Endogenous Mobility)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

            
fp.max_kids_moved  = 90;
fp.min_kids_moved  = 41;
fp.do_counter = 1;
fp.do_counter_type = 2;

pe_skills_moved_end_moving      = NaN(1,fp.max_kids_moved);
pe_skills_receiving_end_moving  = NaN(1,fp.max_kids_moved);

data_baseline = obj_eval.simul_data;

temp_data_baseline  = [];

for j = 1:1:4
    for h = 1:1:fp.n_rip_sim
    temp_data_baseline = [ temp_data_baseline ; data_baseline{j}(:,:,h)];
    end
end

fp.deciles_initial_skills = prctile( log(temp_data_baseline(  temp_data_baseline(:, fp.ind_data.child_age)==9 ,fp.ind_data.child_C)) , [0:10:100]);

set_children_to_move = NaN(length(fp.min_kids_moved:1:fp.max_kids_moved),fp.n_rip_sim);

for h = 1 : fp.n_rip_sim

set_children_to_move(:,h) = randsample(fp.min_kids_moved:1:fp.max_kids_moved,length(fp.min_kids_moved:1:fp.max_kids_moved));

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%draw shocks for receiving families about the probability of 
%leaving the neighborhood as children are moved in;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;

fp.shocks_leaving_neighborhood = rand(fp.schools_distrib(4,3),fp.n_rip_sim);

%Compute the policy effects;
%Loop over the number of children moved in;
ind = 1;
for n_kids =  [1, 5:5:50]
    
    n_kids
    
    fp.n_kids_moved    = n_kids;

    fp.shock_child_moved = set_children_to_move(1:n_kids,:);

    cd(fp.matlab)
    obj_counter2 = obj_function( param_est, fp );
    
    data_counter2 = obj_counter2.simul_data;
    fp.id_moved_children = obj_counter2.id_moved_children;
    fp.id_flights        = obj_counter2.id_flights;

    pe = policy_effects_counter2_end_moving(data_counter2,data_baseline,fp);
    
    pe_skills_moved_end_moving(ind)=pe.pe_skills_moved;
    pe_skills_receiving_end_moving(ind) = pe.pe_skills_receiving;
           
    ind = ind + 1;
end

%Load the original policy effects without endogenous moving decision of
%parents (computed in main/master_counterfactuals);
load('pe_skills_receiving')
load('pe_skills_moved')

grid_x= [1, 5:5:50];
grid_y = grid_x;
cd(fp.paper)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%Appendix Figure E-5: Policy and Scaling Effects on Skills (with Endogenous Mobility)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;

h1=figure;

scatter(grid_x ,pe_skills_moved(1:length(grid_y)) ,  'b')
hold on
p = polyfit(grid_x ,pe_skills_moved(1:length(grid_y)),2);
y = polyval(p, grid_y );
plot( grid_x , y , '--b','LineWidth',2,'MarkerSize',10)
    
scatter(grid_x ,pe_skills_moved_end_moving(1:length(grid_y)) ,  'r')
p = polyfit(grid_x ,pe_skills_moved_end_moving(1:length(grid_y)),2);
y = polyval(p, grid_y );
plot( grid_x , y , '--r','LineWidth',2,'MarkerSize',10)    
    
legend('Policy Effects (No Endogenous Moving)','Quadratic Approx', 'Policy Effects (With Endogenous Moving)' , 'Quadratic Approx' )
xlabel('N. Children Moved');
ylabel('Mean Effect on Log-Skills (Moved Children)');
set(gca,'XTick',0:5:length(fp.min_kids_moved:1:fp.max_kids_moved))
set(gca,'YTick',0.00:0.02:0.24)
set(gca,'YLim', [0.00 0.24 ])
%hgexport(h1, 'Count1_robust.png',hgexport('factorystyle'), 'Format', 'png');
%hgexport(h1, 'Count1_robust.eps',hgexport('factorystyle'), 'Format', 'eps');
    


h2=figure;
scatter(grid_x ,pe_skills_receiving(1:length(grid_y)) ,  'b')
hold on
p = polyfit(grid_x ,pe_skills_receiving(1:length(grid_y)),2);
y = polyval(p, grid_y );
plot( grid_x , y , '--b','LineWidth',2,'MarkerSize',10)
    
scatter(grid_x ,pe_skills_receiving_end_moving(1:length(grid_y)) ,  'r')
p = polyfit(grid_x ,pe_skills_receiving_end_moving(1:length(grid_y)),2);
y = polyval(p, grid_y );
plot( grid_x , y , '--r','LineWidth',2,'MarkerSize',10)    
    
legend('Policy Effects (No Endogenous Moving)','Quadratic Approx', 'Policy Effects (With Endogenous Moving)' , 'Quadratic Approx' , 'Location', 'Southwest')
xlabel('N. Children Moved');
ylabel('Mean Effect on Log-Skills (Receiving Children)');
%hgexport(h2, 'Count4_robust.png',hgexport('factorystyle'), 'Format', 'png');    
%hgexport(h2, 'Count4_robust.eps',hgexport('factorystyle'), 'Format', 'eps');    


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%Figure E-4: Residential Elasticities w.r.t. Mean School Peer Achievement
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;

h3=figure;
plot(1:1:10,abs(elasticity_moving_by_kid_skills),'--bo','LineWidth',3.5,...
                       'MarkerEdgeColor','b',...
                       'MarkerFaceColor','b',...
                       'MarkerSize',7.5)
xlabel('Deciles of Children`s Skills');
ylabel('Prob. Family leaves if Peer Achievement drops by 1%');
cd(fp.paper)    
%hgexport(h3, 'Elasticities_moving.png',hgexport('factorystyle'), 'Format', 'png');
%hgexport(h3, 'Elasticities_moving.eps',hgexport('factorystyle'), 'Format', 'eps'); 
cd(fp.matlab)



cd(fp.matlab)

